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We explore, for the first time, the effect of double fluctuations on both the diagonal and off- 
diagonal self-energy. We use the T-Matrix equations below T c , developed recently by the Zurich 
group (M.H. Pedersen et al) for the local pair attraction Hamiltonian. Here, we include as well the 
effect of fluctuations on the order parameter (beyond the BCS solution) up to second order in U/t. 
This is equivalent to approximating the effective interaction by U in the off-diagonal self-energy. 
For U/t = —6.0, T/t — 0.05, \i/t = —5.5 and A/t — 1.5, we find four peaks both for the diagonal, 
A(n(n /16, 7r/16), uj), and off-diagonal, B(n(n/16, ir/16), uj), spectral functions. These peaks are not 
symmetric in pairs as previously found. In addition: (a) in A(n(n 7r/16), u), the far left peak 
has a vanishing small weight; (b) in B(n(ir/16, 7r/16), ui) the far left and far right peaks have very 
small weights. The physical picture is, then, that the pair physics in the normal phase (T > T c ) is 
still valid below T c . However, the condensation of the e-h pairs produces an additional gap around 
the chemical potential as in BCS, in other words, superconductivity opens a gap in the lower branch 
of a Hubbard-type-1 solution. 

Pacs numbers: 74.20.-Fg, 74.10.-z, 74.60.-w, 74.72.-h 



I. INTRODUCTION 



High-temperature superconductors display a wide range of behaviour atypical of the standard band-theory of 
metals. In the superconducting state these materials become extreme type II superconductors, with a short coherence 
volume, which one might take as an indication of bound pairs. These features indicate that correlation effects 
might be important in understanding the physical nature of these materials. One of the simplest models featuring 
superconductivity and allowing a systematic study of the effect of electron correlations is the attractive Hubbard 
model. Although this model is unlikely to be a microscopic model of high-temperature superconductors, it is likely 
that understanding it will provide insights into the effect of correlations on measurable properties. We adopt the local 
pairing potential in 2-D as a starting point to study the effect of correlations (beyond a simple BCS approach) due 
to the fact that the HTSC materials do not seem to be described by a mean field approach, as it was recognized 
by Randeria et al Furthermore, Puchkov, Basov and Timusk || pointed out that in the HTSC, the Fermi 
surface is estimated to be Ep = 1 — 2 ev, which is not much larger than the energies probed in infrared experiments 
(4 — 300 meV). Such a low Ep could be a reason for violation of a quasi-particle description, i.e., we do not have 
well-defined elementary excitations. This fact may require the use of the local potential with correlations, and indeed, 
the damping obtained from numerical calculations is comparable to the Fermi energy, for small electron concentration. 

In a previous work the effect of electron correlations on some normal-state properties of the attractive Hubbard 
model was studied using a selfconsistent T-matrix formalism, going beyond simple mean-field treatments Q. It was 
found that for intermediate coupling strengths the attractive interaction gives rise to large momentum bound states 
with energies below the two-particle continuum and with a pronounced effect on the spectral properties. Namely a 
splitting of the free band into two, one of which is associated with virtual bound states. Furthermore, a bending in 
the static spin-susceptibility was observed for temperatures just above the phase-transtion 

In a following paper H, the group of Riischlikon used the functional derivative formalism @]. They derived the 
T-Matrix equations in the superconducting phase valid up to second order in the off-diagonal one-particle Green 
functions. The main feature of that communication was that the T-Matrix, appeared both in the diagonal (first order 
in the T-Matrix) and off-diagonal (order parameter plus a second order contribution in the T-Matrix) part of the 



self-energy. The authors of Ref. || implemented a low order approximation with full fluctuations in the diagonal part 
of the self-energy while using a mean field approximation for the off-diagonal part of the self-energy, i.e., equal to the 
order parameter. 

In this communication, we implement a higher order approximation to the off-diagonal self-energy by allowing 
fluctuations in the order parameter also, i.e., making T(x, x') » U S(x — x'). We keep the full T-Matrix in the diagonal 
part of the self-energy as in Ref. Q. The goal of the present work is to study the stability of the physical picture 
given in Ref. [^), which can be summarized as the appareance of four peaks, symmetric in pairs, both in A{k, k,u)) and 
B(k,k,uj). We caution the reader that we are exploring the effect of double fluctuations by approximating T(xx') 
by US(x — x') in £12 (x, x') as a step forward in taking into account the full effect of correlations. In Section y we 
present the model and the results of Ref. Q . 

In Section III we keep the low order approximation taken in Ref. Q and implement our second order approximation 
in the off-diagonal self-energy. Our implementation is performed with the fast-Fourier transform (FFT) and we 
disscuss the physical meaning of our results. Section IV concludes. 



II. THE MODEL AND THE T-MATRIX EQUATIONS. 

The Hubbard Hamiltonian is defined as 

ff=-* Xi c \a c v° + u ^2 n iT n H ~ V^m* , (1) 

<W>a I la 

where the c\ a {ci„) are creation (annihilation) operators for electrons with spin a. The number-operator is ni„ = c\ a ci ai 
t is the hopping matrix element between nearest neighbours I and I', U is the onsite interaction and /i is the chemical 
potential in the grand canonical ensemble. Here we consider an attractive interaction, U < 0. For a review of the 
attractive Hubbard model see Micnas et al 0. Previous authors have used this model to study the bismuthate 
superconductors 

By starting with the Nambu Green function 

Gxx'(t,t') = -(t t [*x(r) (8*1(^)1) , (2) 



where ® means the tensor product and T T means the time ordering of the two Nambu operators, where these are 
defined, at position x and imaginary time, r, by 



CxT 

J 



the authors of Ref. || write the solution to Dyson's equation in k, w-space as 

(4) 



Gn(k,iuj n ) = -G 22 (k, -iu n ) = 
iu) n + £ k - S 2 2(k,iw„) 



(ioj n - e k - Sn(k,iw n ))(jw„ +e k - S 2 2(k,«w T ,)) - Si 2 (k, w n )E 2 i(k, iu n ) 

Gi 2 (k,iu n ) = C?2i(k, -iw n ) = 

Si 2 (k,io;») 

(iui n - Ek - £11 (k, iuj n ))(iuj n +£k - S 22 (k, iu n )) - T, 12 (k, icj„)E 2 i(k, iw n ) 



(5) 



where uj n = 7r(2n + l)//3 are the fermionic Matsubara frequencies and (3 = l/(ksT) the inverse temperature. The 
d-dimensional dispersion is given by £k = — 2t^ a=1 cos(fc Q a Q ) — fi. S is the self-energy matrix Q and d is the lattice 
dimension. 

Then, by using the functional derivative techique Q , the authors of Ref. (5| get the self-energies to second order in 
Gi 2 (G 2 i) 

E u (x, x') = G 22 (x, x')T(x', x) + G 21 (x, a)G 22 (a, b)T{b, x)Gi 2 (6, x')T(x', x) , 
Ei 2 (x, x') = A(x)tf (x - x') + G 22 (x, a)Gi 2 (a, b)T(b, x)G„ (b, x')T(x', a) , 
£21 (x, x') = A'(x)«J(x - x') + Gn (x, a)G 2 i (a, b)T(b, x)G 22 (b, x')T(x', a) , 

E 22 (x,x') =G 1 i(x,x')T(x',x) + Gi 2 (x,a)G 1 i(a,6)T(&,x)G 2 i(6,x')T(x',x) . (6) 
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where T(x, y) in reciprocal space has the following form 



with 



T(q,ie m ) = - — T — : — r , (7) 

1 - Ux(q,i£ m ) 



x(g) i£ m ) = j-fg ^ G 22 (fc - q, iuJ n - fem)G'ii(fe) iu n ) , (8) 



and e n = 27m//? are the bosonic Matsubara frequencies. This solution (see Eq. (H)) is valid for A/W <C 1, where 
W = 2cft is the bandwidth. In Eq. (^|) summation over repeated indices is understood (space and imaginary time). 
The physical meaning of Eq. (P) is the repeated scattering of two particles many times, without polarization of the 
medium. This is equivalent to the T— Matrix approximation which is valid for small carrier concentrations. Another 
approximation involved in the derivation of Eq. @ is the condition that A(T) << W . This is equivalent to a 
perturbation expansion to second order in A(T)/W. For a solution of the Eqs. ([|,[7]|8|(]) one would also need to fix 
the chemical potential from the particle number using 

P(T,H)= lim 52 G(k,iu) n )exp(iu n ri), (9) 

where p is the electron concentration per spin and is defined in the interval [0, 1]. Thus, the set of Eqs. (||,0,|] |]) 
represents a set of non-linear self-consistent equations which needs to be solved numerically. 

We note that an expansion of the final equations, Eqs. @, to first order in U simply gives the wellknown BCS 
expressions. To second in U the result is identical to that of Martin- Rodero and Flores plfl . 

III. NUMERICAL RESULTS 

Before we disclose our approximation and present a numerical solution to it, we will discuss the approximation 
made in Ref. Q. The authors of this reference approximate the set of Eqs. (||) in the following form |l2|,13| 

£n(x,x') = G 22 (x,x')T(x',x) 
S 12 (x,x') = A(x)5(x-x') 
S 21 (x,x') = A*(x)5(x-x') 

E 22 (x,x') = Gii(x,x')T(x',x) , (10) 

which is equivalent to keeping full fluctuations in the diagonal part of the self-energy while performing a mean field 
approximation in the off-diagonal self-energy, Si 2 (x, x') = A(x) S(x — x'). Remember fluctuations enter through the 
T-Matrix, and this has been set equal to zero in the off-diagonal self-energy. The idea behind the point of view taken 
by these authors is that the effect of correlations is mainly present in the diagonal self-energy, i.e., correlations are 
linear in the T-Matrix. In addition, the authors of Ref. Q took the pragmatic view of taking one step at a time, i.e., 
to study correlations in the same manner as they were studied in the normal state. 

Here, we want to include fluctuations in the order parameter, or equivalently in the off-diagonal part of the self- 
energy. For this, we make T(x, x') « U5{x. — x') in the off-diagonal self-energy while keeping the full fluctuations in 
the diagonal self-energy. Then, our approximation is equivalent to the following set of equations 

E 11 (x,x') = G 22 (x ) x')T(x',x) 
Ei 2 (x, x') = A(x)(5(x - x') + f/ 2 Gn(x, x')Gi 2 (x', x)G 22 (x, x') 
E 21 (x, x') = A*(x)(5(x - x') + f/ 2 G 22 (x, x')G 21 (x', x)Gn(x, x') 

S 22 (x,x') = G 11 (x 1 x')T(x',x) (11) 

which are local in real space where they can be easily evaluated numerically. The Greens functions are determined 
selfconsistently using Dyson's equation where the order parameter is used as input (see below). The technical aspects 
of the numerical solution using the FFT-technique has been detailed in S. We consider that the set of Eqs. ( p"l| ) 
contains double fluctuations since we have full fluctuations in the diagonal part of the self-energy and fluctuations 
in the off-diagonal self-energy. The latter come through the one-particle Green functions themselves which must be 
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calculated self-consistently. The reason that we pursue the present approximation is that we believe that including 
fluctuations in the off-diagonal self-energy will have a strong influence on the dynamical properties. The outcome of 
our numerical calculations vindicate this. Although the choice of T(x, x') w US(x — x') in the off-diagonal self-energy 
does not fully take into account all the fluctuations, it keeps the numerical implementation under control with present 
time Workstation facilities. We leave to the future an implementation of the full T-Matrix contribution. 

The numerical simulations were performed in two dimensions (2D), for U/t = —6.0, fi = —5.5, A = 1.5 and 
T/t = 0.05. We have fixed the order parameter and the chemical potential because computationally is much easier 
to find one root rather than two. The selected value of A/t = 1.5 for U/t = —6.0 is inside the range of validity of the 
approximation and close to the numerical value found in Ref. Q . As is explained in Ref . |l4| numerical convergence is 
based on the fact that the diagonal self-energy depends neither on [i nor A. This assumption is not longer valid here, 
making the program more complex. Because of these reasons, we have decided to fix A and [i and find the carrier 
concentration, p. In addition, there is more than one solution for the equation p(A, [i) = const (see Fig. 8). Here one 
must note that due to the Mermin- Wagner theorem no phase-transition is expected in 2D systems with a continous 
symmetry, and the formalism employed is therefore too simple to describe a Kosterlitz-Thouless phase-transition. 
Nonetheless, it is observed that the formalism does give a phase-transition, in agreement with the results above the 
transition temperature, where the signaling of a divergence of the T-matrix was observed. 

In Figure 1, we present the diagonal one-particle spectral function, A(k,Lo), defined by 

A(k,w) = — lim Im[G n (k,w+i5)l (12) 

The value of the gap, A/t = 1.5 is within the region of validity of the expansion, Eqs. (^), i.e., A/W < 1, with 
W = 8t the bandwidth in 2D. We observe three visible peaks, two symmetric around the chemical potential, \i, and 
the third is the the upper Hubbard branch. However, there is a fourth peak for lo < 0, which for n = is around 
lo k, —3.8. The equivalent peak for lo > is around ~ +3.25. The peak atw« —3.8 has a vanishingly small weight 
(see the inset figure). So, we conclude that the symmetric peak structure does not hold anymore (at least for the 
most extreme peaks). If we consider that the peak at lo w —3.8 is indeed small then we are left with effectively three 
peaks. These results are different from the ones obtained by the Zurich group Then, the picture which emerges 
from A(h, lo) is the appearence of three energy branches, a higher one which is already present in the normal phase 
and two lower ones which are symmetric around the chemical potential. These two superconducting lower branches 
correspond to the lower energy branch of the normal state. Now superconductivity (A ^ 0) opens up a gap in the 
normal lower energy branch. These two superconducting low branches are, then, due to the pairing of the electron 
- hole pairs accross the chemical potential, similar to the BCS case. The only difference with respect to the weak 
coupling (BCS) limit is the appearence of the higher energy branch. 

In Figure 2, we present the off-diagonal one-particle spectral function, B(k, lo), which is defined as, 

B(k,u>) = — lim Im[G 12 {k,u; + i5)], (13) 

From Figure 2 we see that this function has two symmetric and visible peaks. In the inset we show a blow up of 
the energy scale. We observe two additional two peaks at lo w +3.9, —4.1, for n ~ 0, which are not symmetric and 
their weights are small. So that we argue that they can be neglected. Finally, there are two symmetric peaks around 
fj, (lo = 0). The conclusion to be drawn from these observations is that the off-diagonal spectral function remains 
qualitatively equal to the the BCS case, i.e., the inclusion of fluctuations both in the diagonal and off-diagonal 
components of the Nambu self-energy does not change the basic BCS condensation picture. Only the diagonal 
spectral function suffers the effect of correlations (Figure 1), as should be, due to the presence of strong correlations. 
We have chosen k x and k y in units of 7r/16. 

In Figure 3 we present the imaginary part of the two-particle Green's function, — 7roG2(q, lo). The Cooper resonance 
almost dissapears since the peak weights for lo < are negligible (see inset). In addition, for every value of k we 
find (lo > 0) a main peak with some two broad peaks of reasonable weights. These additional broad peaks have an 
effect on S(k, iu> n ). As S(k, iLo n ), G(k, iLo n ) and F(k,iu> n ) are coupled together, then, these peaks which in the are 
going to have an effect on the quasi-particle spectra. All this is due to the fact that our equations are self-consistent. 
We observe a well defined two-particle structure for all values of momentum for lo > 0. The fact that the Cooper 
resonance is, in certain sense, washed out is probably due to the fact that the effect of double fluctuations favor 
pairing formation, which implies that the interaction is stronger. A discussion of this point is treated in the paper of 
Ranninger and Robin p5|. Another difference between our results and those of the Zurich group, is that our band of 
two particles is wider than theirs. This has deep physical consequences, as we will see later. 

In Figure 4 we present the imaginary part of the diagonal self-energy, — Im[T,(n(ir/16, 7r/16), lo)] for the same 
parameters of Figure 3. We observe an opening of the gap at the chemical potential, lo = 0. This gap is much bigger 
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than twice the order parameter. In the inset we have a blow up of the energy scale concluding that it is composed 
of two peaks, one of them broad. The authors of Ref. ||]4| found that ImS almost has a peak in the normal phase. 
The presence of this peak is the origin to the two peaks in A(n(7r/16, 7r/16), uj)]. In the superconducting phase, the 
authors of Ref. j^] find that ZmS still has one peak which implies the presence of four peaks, symmetric in pairs, 
both in A(n(ir/16, 7r/16), uj)] and B(n(7r/16, 7r/16),u;)]. Now, in our approach with double fluctuations, we observe 
that ImE has no well defined single peak structure. This, we argue, will wash out the pair symmetry. The presence 
of double fluctuations is responsible for the cancellation of the leftcst peak of A(n(Tr/ld>, tt/16), uj)]. So, the picture 
which emerges is a Hubbard-type-I solution combined with a superconducting gap in the lower branch ^M . 

In Figure 5 we show the real part of the diagonal self-energy, i?e[E(n(7r/16, 7r/16), uj)], for the same parameters 
given in Figure 3. We observe in Figure 4 that the width of these peaks is bigger than those of the Zurich group. The 
consequence of this is that we are no longer able to define a well define peak in the self-energy, so that the four peak 
structure (with symmetric pair of peaks) of the Zurich group is not longer valid. Thus, the picture is the following: 
above T c we have two peaks which are due to the pair formation physics, i.e., the upper and lower branches; below T c 
one of the bands is split into two to the spontaneous symmetry breaking, A ^ 0. 

In Figures 6 and 7 we show the frequency dependence of the real and imaginary part of the off-diagonal self-energy, 
i.e., i?e[Ei2(?i.(7r/16, tt/16), uj)] and — /m[Ei2(n(7r/16, tt/16), uj)], for the same parameters as given previously. For the 
case of — 7m[Ei2(n(7r/16, 7r/16), uj)] we have not shown the full energy range. However, it is fully antisymmetric. 
To the best of our knowledge, this is the first presentation of these results in the literature. We observe that 
7m[Ei2(n(7r/16, 7r/16), uj)] and Im[T,(n(n /16, tt/16), uj)] have almost the same values for large frequencies. However, 
the position of their peaks is at a different frequency for each n. In addition, the effect of correlations (fluctuations) 
is important since they decrease the value of i?e[Ei2(n(7r/16, it/ 16), uj)] from A down to almost A/2, for n = and 
uj = 0. If we take a closer look at Eqs.([|]|) we see that ^12 plays the role of an energy gap. Then, the energy gap 
is (k, w)-dependent and it of the same order of magnitude than the order parameter, A(T). This analogy leads us 
to conclude that the order parameter, A(T), and the energy gap are two different quantities. So, a local interaction 
produces a (k, w)-dependent energy gap, which is the quantity that experimentalists are most likely measuring. 

In Figure 8 we have a plot of p vs p for different fixed values of the order parameter, A. We see that there is an 
abrupt change of density as a function of p. Also, for a fixed value of p there are up to three different values of the 
chemical potential. Then, our scheme is not valid for large densities, a condition we take as granted by the definition of 
the T-Matrix formalism. On the other hand, when we increase the temperature the step of p vs p tends to dissapear. 
This implies that our scheme is valid for temperatures which are not too low. We would like to point out that the 
implementation of the T-Matrix approach in the normal phase (A = 0) did not converge for small temperatures and 
high densities ||. We interpret the non-convergence of the T— Matrix program above T c as due the fact that we are 
approaching a parameter region where there are no physical solutions. For example, in another context, Figueira and 
Foglio p7| , have found multivalues for the chemical potential in the periodic Anderson model and they choose the 
solution of minimun Helmholtz free energy. 



IV. CONCLUSIONS. 



Using the analytical results of Ref. |3| , we have implemented a program where we have included fluctuations both 
in the diagonal and off-diagonal self-energy. This approach we call double fluctuation since we have full fluctuations 
in the diagonal part of the self-energy while going beyond the mean field solution (BC'S) in the order parameter. 
This represents a higher order approximation than those used in Ref. We have presented the one-particle and 
two-particle spectral functions pointing out the differences with the low order approximation performed by the Zurich 
group (5j . We have indicated that the pair formation, which is a phenomenon valid even above T c remains valid in our 
scheme. However, the effect of fluctuations in the off-diagonal part of the diagonal self-energy (double fluctuations) 
does not produce a well defined peak in the diagonal part of the one-particle self-energy and this breaks the four peak 
structure found by the Zurich group, resulting in a clearer picture since only one of the Hubbard bands is split into 
two around the chemical potential. The off-diagonal spectral density is qualitatively similar to the one of BCS, but 
the diagonal spectral function suffers the effect of strong correlations. As has been said before, the physical picture 
which emerges from our work is that of superconductivity opening a gap in the lower branch of a Hubbard-type-I 
solution. We have performed additional calculations and found that the opening of the correlation gap is hard to 
achieve. For U/t = —12.0, T/t = 0.1 and p/t = —8.5 we find that the correlation gap is more or less visible but 
not fully developed. This is due to the fact that the shifting of the electronic states around the chemical potential 
sends electronic states above p and this closes the correlation gap. We have also discussed the parameter space where 
our approximation is valid and we have made a connection with the calculations above the transition temperature 
discussing the validity of each. To end, we mention that a previous calculation by Rodriguez - Nunez, Cordeiro and 
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Dclfino |n| has used the ideas presented in this work within the framework of the moment approach of Nolting |T^] 
in the superconducting phase. 
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FIG. 1. The diagonal one-particle spectral function, A(n(pi,pi), lo) vs lo for different momenta along the diagonal of the 
Brillouin zone (k = (n,n)n/16) for U/t = —6.0, T/t = 0.05, A = 1.5 and fi = —5.5. We have used an external damping of 
5/t — 0.1. We have used 16 x 16 points in the Brillouin zone and 1024 Matsubara frequencies. After self-consistent calculation 
of the coupled non-linear equations, we obtain a density, p = 0.035942. we have runned our source code in single precision 
requiring 18 MB of RAM memory. Each iteration takes 2.5 minutes of CPU time. 

FIG. 2. The off-diagonal one-particle spectral function, B(n(7r/16, 7r/16), w) vs lo for different momenta along the diagonal 
of the Brillouin zone. Same parameters as in Figure 1. 

FIG. 3. — I m[G 2 {m(ir/ 16, n/ 16), lo)] vs lo for different momenta along the diagonal of the Brillouin zone (q = (m, m)7r/16). 
Same parameters as in Figure 1. 

FIG. 4. — 7m[E(n(7r/16,7r/16),a>)] vs lo for different momenta along the diagonal of the Brillouin zone (k = (n,n)n/16). 
Same parameters as in Figure 1. 

FIG. 5. .Re[E(n(7r/16, 7r/16), lo)] vs lo for different momenta along the diagonal of the Brillouin zone (k = (n, n)7r/16). Same 
parameters as in Figure 1. 

FIG. 6. — Jm[Ei2(n(7r/16, 7r/16), lo)] vs lo for different momenta along the diagonal of the Brillouin zone (k = (n, n)7r/16). 
Same parameters as in Figure 1. 

FIG. 7. iie[Ei2(n(7r/16, 7r/16),a>)] vs lo for different momenta along the diagonal of the Brillouin zone (k = (n, n)7r/16). 
Same parameters as in Figure 1. 



FIG. 8. p vs n for different fixed values of the order parameter, A, at T/t = 0.05 and T/t = 0.2. 
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